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ABSTRACT 

Often, in solving an elliptic equation with Neumann boundary conditions, a compatibility 
condition has to be imposed for well-posedness. This condition involves integrals of the 
forcing function. 

When pseudospectral Chebyshev methods are used to discretize the partial differential 
equation, these integrals have to be approximated by an appropriate quadrature formula. 
The Gauss-Chebyshev (or any variant of it, like the Gauss-Lobatto) formula can not be used 
here since the integrals under consideration do not include the weight function. A natural 
candidate to be used in approximating the integrals is the Clenshaw- Curtis formula, however 
we show in this paper that this is the wrong choice and it may lead to divergence if time 
dependent methods are used to march the solution to steady state. 

We develop, in this paper, the correct quadrature formula for these problems. This 
formula takes into account the degree of the polynomials involved. We §how that this 
formula leads to a well conditioned Chebyshev approximation to the differential equations 
and that the compatibility condition is automatically satisfied. 
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INTRODUCTION 

We deal here with a problem encountered in the solution via Chebyshev spectral colloca- 
tion discretization of the Poisson equation with homogeneous Neumann boundary conditions. 
The problem arose in the context of solution of the pressure Poisson equation in a time-split 
algorithm for the incompressible Navier-Stokes equations. For this problem to be well-posed, 
the source term must satisfy a compatibility condition - , the numerical analog of this condi- 
tion using straightforward Clenshaw- Curtis quadrature formulae was found to be numerically 
ill-conditioned and produced large distortions in the spectral solution. 

In Section I, we outline the time-split algorithm and describe the difficulty encountered; 
we proceed in Sections II and III to analyze the difficulty, first in terms of an equivalent 
parabolic equation which would be utilized, for instance, in the iterative solution of the 
Poisson equation of interest, then in terms of the steady equation itself. The proper quadra- 
ture formula is also developed, which alleviates the numerical difficulty. In Section IV we 
show numerical examples to demonstrate the numerical inconsistency and its resolution. 

I. TIME-SPLIT ALGORITHM 

A splitting method is employed in many simulations to advance the solution of the incom- 
pressible Navier-Stokes equations from time t n to f" +1 . Writing the Navier-Stokes equations 
in vector notation, with u representing the velocity (u,v,w) we have: 

u ( + u ■ Vu = - VP + vV 2 u (1.1a) 


and 


V • u = 0 


(1.16) 


in the domain fi, and : 

u = 0 on the boundary T 

The split scheme first advances u n to an intermediate solution u* by solving: 

Uj + u* • Vu* = i/V 2 u* 


( 1 . 2 ) 


u* = g* on T 


The intermediate boundary condition u* = g* is discussed in [1]. Finally, the solution is 
advanced from u* to u n+1 via: 

u” +1 = -VP" +1 

V • u" +1 = 0 (1-3) 

h ■ u n+1 = 0 on T 
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where h is the unit normal to the boundary I\ Note that the final, “pressure correction” 
step by itself is a set of inviscid equations; and is well-posed when boundary conditions on 
the normal component of the velocity only are enforced. At the end of the full step there 
exists a non zero tangential component of velocity on the boundary. The magnitude of this 
slip velocity ~an be reduced to 0(At 3 ) by a proper choice of the intermediate boundary 
condition on u* [1]. 

Using backward Euler time discretization for Eq. (1.3) yields: 

u "+i _ u * _ AiVP n+1 (1.4) 

The pressure step is actually carried out in two parts. First, the divergence of Eq. (1.4) 
yields: 

V 3 P n+1 = _Lv-u* (1.5) 

where V • u n+1 — 0 is enforced. Then the velocities are updated using Eq. (1.4). 

Note that this formulation requires a boundary condition for the pressure. This poses a 
problem, since there is no natural boundary condition for pressure. The use of a condition 
derived from enforcing the normal momentum equation at the boundary was attempted, but 
was apparently inconsistent as it resulted in explosive instability. Fortunately, analysis of 
the split scheme itself yields a self-consistent pressure condition: 

n • VP n+1 |r= 0 (1.6) 

since both n-u* and n- u" +1 are zero on the boundary. The error involved in this specification 
is, we believe, related to the overall splitting error of the scheme. It is also known that the 
error due to imposition of an inconsistent Neumann pressure boundary condition is isolated 
to a thin “boundary layer” [4] . 

It is easily shown using the divergence theorem that 

J (V • u*) dft = 0 (1.7) 

is required for the pressure Poisson equation to be well-posed with this boundary condition. 
For application of this algorithm in closed-system problems with homogeneous Dirichlet 
velocity boundary conditions, numerical tests indicate Eq. (1.7) is satisfied in general. How- 
ever, in applications of a recently- developed non- reflecting outflow boundary treatment [5] 
to simulations of flow-through systems, it was found that Eq. (1.7) did not hold discretely 
for the intermediate velocity field. As will be shown in the numerical examples of Section 
IV, this leads to large distortions in the computed pressure field. It was decided that since 
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there is no solution to the problem: 


V 3 P = cast ^0 in fl 
VP • ra = 0 on T 

the pressure Poisson equation would be modified to read: 

Vu*-^(V. u *)dfi] (1.8) 

The integral in Eq. (1.8) was implemented using the straightforward Clenshaw-Curtis 
quadrature formulae, appropriate for Chebyshev collocation. Although the pressure field 
distortions were much reduced, their magnitude was still unacceptable. The analysis pre- 
sented in the following show the inconsistency in the use of these formulae, and the proper 
quadrature formula to recover a smooth solution to Eq. (1.8). 
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II. TIME DEPENDENT PROBLEMS 

Consider the parabolic equation 


u t =u xx +g 

(2.1u) 

u(x, 0) = 0 

(2.1*0 

with the Neumann boundary conditions 


«*(M) = «*(-!,<) = 

(2.2) 


We assume also that 






j g(x)dx = 0 

(2.3) 

so that 

1 




jtL u{x - 

,t)dx = «j.(l) — u x (— 1) + J ^gdx = 0 

(2.4) 

and therefore in view 

of (2.1b) 





J u(x, t)dx — 0 for all t. 

(2.5) 


In the pseudospectral Chebyshev method we seek an x polynomial u^(x,t), that satisfies 
the boundary conditions (2.2), such that (2.1a) is satisfied at the points Xj = cos namely 
we seek u N(x,t) such that 


duN 

~di 


d 2 Uff t 4 _ 

-faT + 9 at x ~ X 3 




uw(x,0) = 0 


(2.6a) 

( 2 . 66 ) 


and 


£(-U)-£<M)-o. 


(2.7) 


We refer the reader to [3] for a discussion of the stability and convergence of (2.6) and (2.7) 
to (2.1) and (2.2). 


Here we are interested in the question whether the numerical approximation satisfies 
the conservation property (2.5). Since the numerical solution u#(x,t) defined in (2.6) is a 
polynomial of degree N in x it is natural to consider the Clenshaw Curtis quadrature formula 
[2]. It uses the collocation points of (2.6a) and it is exact for polynomials of degree N. 
Lemma (1.1). Let N be an even number and let qj be defined by 


a. 


__ 4 * 1 


2 irlk 1 . . r 

cos • — 0 < l < N 

N c fc 


( 2 . 8 ) 
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a 0 = a N = 


N 2 - 1 

Then for any polynomial h(x ) of degree at most N 


where c 0 = cjv = 2 c* = 1 0 ^ k ^ N 


/: 


1 W 

h(x)dx = ^2 h(xi)on 

1=0 


(2.9) 


where ®j = cos 

A close inspection of (2.6a) reveals that, in general, ujv (®,<) does not satisfy the con- 
servation condition (2.5). Indeed, since g(x) in (2.1) is not a polynomial of degree N, the 
quadrature formula (2.9) is not exact for g, whereas it is exact for the other terms in the 
equation. 

To remedy the situation it is customary to modify (2.6a) by seeking Uff(x,t) that satisfies 


6un 3 2 Un 1 ^ , ^ A 

~W + ~dx 2 ~ = 9 ~ o^ 9 ( Xl ' ai at x = x i 1 < 3 < N 


(2.6c) 


/=0 


where qj are defined in (2.8). Of course uj v(x,t) still satisfies (2.1b) and (2.2). Thus the 
right hand side of (2.6c) has zero mean if the Clenshaw- Curtis formula is being used. Unfor- 
tunately, even with the above modification the solution of (2.6c) does not satisfy the discrete 
analog of (2.5); in fact we can state: 

Lemma (1.2): Let u^(®,i) be the solution of (2.6c,b) with the Neumann boundary condi- 
tion (2.7). Let gw(t) be the coefficient of the JV’th Chebyshev polynomial in the expansion 
of g(x, t), namely 

9n 


Then 


where 


un{ 


1 g{x)T N {x) J _ 

(2.10) 

;) = cos AT(cos -1 ®). 

(2.11) 

: N 2 _ i {MO l SN{t)dr } 

(2.12) 

i r 1 ^(®,t)iv(x)^ 

7T 7-1 y/l — x 2 

(2.13) 


Remark: Equation (2.12) points out the difficulty in the method (2.6c). Suppose for exam- 
ple that g is a time independent function, then the solution of (2.1a,b) converges to a steady 
state solution and therefore its highest coefficient itw(Z) is a bounded function of t. However 
the second term in (2.12) is not; in fact 

f 9n(t )dr = tg N 
Jo 
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and this diverges linearily in t. 

The above does not contradict the usual stability argument since for fixed t, 
converges to tt(x,i) as N tends to infinity , but the usefulness of (2.6c) as a method for 
marching to steady state is doubtful. This will be demonstrated further in the numerical 
experiments. 

Proof: We denote by PNg(x,t) the (unique) interpolation polynomial of degree N in x that 
interpolates g(x , t ) at the points x, = cos ~ 0 < l < N. With this notation equations 

(2.6c)-(2.7 ) can be rewritten as 

IT = + PnS " + (j4 * + b)t ' n(x) (2 ' 14o) 


u w (i,o) = o = (M) = 


(2.146) 


By comparing the coefficients of the highest Chebyshev polynomial in both sides of (2.14a) 
one gets 


A = 


ditw 


dt 


— 9n 


N' 


(2.15) 


Multiplying now (2.14a) by otj and summing we get 


4 X> K (X„ i)ay = £ fiat + £ g(xt)at - £ + £(.4*; + B)T' N (xj) 

dt 1= o 1=0 j= 0 J=0 J-0 

(2.16) 

We use now the exactness of the Clenshaw Curtis formula and the fact that Tjj (x^) = 0 1 ^ 

j < N — 1 to get 


K (xt,t)at = - ^(- 1.0 + UN> ao 


2N 2 
N 2 - 1 


dittf 

~dt 


- 9n{1) 


N 


(2.17) 


Integration of (2.17) yields (2.12) and the proof is completed. 

It is obvious that the problem lies in the use of the Clenshaw Curtis quadrature formula. 
Specifically, we need to use an accurate quadrature formula that does not use the boundary 
points. Such a formula is constructed in Lemma(1.3). 

Lemma (1,3): Let /(x) be a polynomial of degree N — 2 at most. Let Xj = cos ^ then 


/ f(x)dx = 9{xj)P; 

J ~ 1 j = 1 


(2.18a) 


where 
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(2.186) 


Pi = (-l) i+1 - 


+ a j j = l,...,N -1 


N 3 - 1 

and aj are defined in (2.8). 

Proof: Since f(x) is a polynomial of degree N — 2 , it is uniquely determined by its values 
at the interior points Xj,j = l,...,N — l. In fact 


/(*) = E H x i) t '- x ' > 

3 = 1 


X~Xj TnM’ 


Therefore 


where 


/ f(x)dx = Y, f( x j)Pj 
J ~ l j=l 

Pi = f 

J- 1 X — Xj 1 ff{Xj) 


We use now the Clenshaw Curtis formula to evaluate the integral (2.20). Thus 

q _ Y' T n( x i) a i | 

Pi ~ JL, _ -x -T"(x ) + J 

{= o X 3 1 N\ X 3) 

‘*3 

_ 7U- 

- (l-x j )T^) Qo+ (-l-^M(x,0 + ’ 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


= (- 1 ) 


i+i 


W 2 - 1 


+ a 3 


and the quadrature formula (2.18) is thus established. 

Using Lemma (1.3) we suggest the following algorithm for the pseudospectral Chebyshev 
discretization of (2.1)-(2.2). In the new algorithm we seek ujv(x } t) such that 


duff d^ut4 
dt dx 3 


N - 1 


+ J-rE 9( x j)Pj at x = x i 3 = !» • • • > N ~ 1 

1 3=1 


(2.22o) 

(2.226) 

(2.2c) 


it/^(x, 0) = 0. 

In the new method the compatibility condition (2.5) is satisfied. In fact we state 
Theorem (1.1): Let ujv(x,f) be the solution of (2.22). Let (3j be defined in (2.21) then 
u N (xj, t)(3j = 0 for all t. 
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Proof: We multiply (2.22a) by {3, and sum to get 


j N - 1 N - 1 Q2„ N-l N-l 

± Y, m = E Vr (*i. ■ Oft + E *(*;)& - E (2- 23 ) 

df j=i j= i ax i=i i=i 


We use now Lemma (1.3) and the fact that is a polynomial of degree N — 2, to get 

N ~ x d 3 u N , . f l d 3 u N j du N ,. du. 

i'=i 

and therefore 


^ l 0*ujir, jNjtJ /* 5 3 uaTj ^ 0«jv, ^ 

E -ajr* = - -g^(-M) = 0 


(2.24) 


and since u^(x,0) = 0 


which proves the theorem. 


d N_1 

— E u N(Xj t t)Pj = 0 


N - 1 

53 «w(*i» t)0j = 0 

J =1 


(2.25) 
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III. STEADY STATE PROBLEMS 

In this section we discuss the equation 


v-xx — 9 (3.1c) 

u x (±l) = 0 (3.16) 

Equation (3.1) can be viewed as the steady state version of (2.1a), however this time we 
need one more condition to get a unique solution. We impose the condition 

J u dx = 0. (3-2) 

In order for (3.1a) and (3.1b) to be compatible g(x) has to satisfy the condition 

J g(x)dx = 0. (3-3) 

We have demonstrated in Section II that the use of the Clenshaw-Curtis quadrature 
formula (2.8) causes problems in solving (3.1) by trying to evolve the solution of (2.1) to 
steady state. Here we would like to discuss the direct solution of (3.1) via the influence 
matrix technique. We will demonstrate that the quadrature formula (2.18) developed in 
Section II should be used, rather than the Clenshaw-Curtis formula (2.8). 

The solution of (3.1)-(3.2) via the influence matrix technique involves seeking the ap- 
proximation of (3.1) as a sum of the solution of three problems: 

Problem 1: Seek a polynomial u>#(x) such that 


d 2 w 


' N , i / • ✓ AT 

, j- = 9- 9o at x = Xj = cos — 1 < J < N 
w N (\) = 0 = iwjv(-l) 

and go is an approximation to / g{x)dx. We will consider either 

N 

9o = £</(*>; 
j = o 

aj given in (2.8) 


(3.4) 


(3.5c) 


or 


N-l 

9o = £ 9{xj)Pi- 

j=i' 

f}j given in (2.18b) 

Problem 2: Seek a polynomial v r (x) such that 

cPv 1 


(3.56) 


dx 2 


= 0 at x = Xj j = 1, . . . , N — 1 


(3.6c) 
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(3.66) 


t/(-l) = 0 v J (l) = 1 


Problem 3: Seek a polynomial v n (x) such that 


d x 3 


= 0 at x — Xj j = 




(3.7a) 


v n (-l) = 1 t> J/ (l) = 0 


(3.76) 


Given wjftV^v 11 we look for a solution of the form 


un(x ) = u;Ar(x) + diV T (x ) + d 2 v Il (x). 


(3.8) 


Clearly Um satisfies 


d 2 Uff(i c) 


dx 2 


= g(x ) — go at x — Xj 1 < j < N 


for any cfi , cia- The constants d\,d 2 are determined by the boundary conditions (3.1b). 
In our case (3.1b) (3.2) imply 


0 = £<-»> = is P 1 + + <^(-D 

0 = = *»»(!) + + ^^(1) 


(«) 


dx 


dx 


dx 


dx 


0 = J Ufif(x)dx = J Wfj(x)dx + d\ J v\x)dx + d 2 J v n (x)dx 


(6)(3.9) 

(c) 


Ideally (3.9a) and (3.9b) should yield the same equation for di, d 2 . The integrals in (3.9c) 
have, or course, to be discretized. 

In this simple model problem we are able to evaluate explicitly w^(x),v I (x), and v n (x). 
In fact we can state 

Lemma (2.1): The solutions of problems 2 and 3 are given by 

v' = «" = ^ (3.10) 

Proof: v I (x),v Il (x) in (3.10) are the only polynomials that satisfy (3.6), (3.7) respectively. 

In order to get an expression for wm(x) we introduce the interpolation polynomial Pm -29 
that interpolates g(x) at the points Xj 1 < j < N — 1. By (2.19) 


9(x)= Yl9(^j) TN ^ 
j=i 


-Xj Ttf(xj)' 


(3.11) 
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(3.12a) 


It is clear form (3.4) that wn(x) satisfies 

d?w n 


= Pn- 29 ~ 9o 


'U'w(l) = 0 = wpt(— 1). (3.126) 

Note that equation (3.12a) holds for any — 1 < x < 1, not just for the grid points Xj ; it is 
therefore possible to integrate (3.12a) to get 


= ^ (_1) + Ij PK - i3 ^ ) - (313) 

and 

w N (x) = -^p(-l)(l + x ) + J x (® - v)[PN-2g{v) - 9 o]dr]. (3-14) 

Substituting now the condition ttfw(l) = 0 we get 

^r (-1) = ~\ Ij 1 ~ - *]<*>• (3-15) 

We are ready now for the main result of this section. 

Lemma (2.2): Let g 0 be defined by the use of the Clenshaw-Curtis formula as in (3.5a), 
then equation (3.9a) is incompatible with (3.9b). 

Proof: Using the expressions (3.10) for v 1 and v 11 in (3.9a,b) we get 


(3.14) 


(3.15) 


dwff . . _ dj. ^2 

dx { 2 2 


dw n . . _ di d 2 

HhT^ ' ~ 2 ” 2 ' 


(3.16) 


However from (3.13) 


- s°w- 


(3.17) 


Since P N -29 is a polynomial of degree N — 2, the quadrature formula (2.18) is exact and 
therefore N 

[ [Pn- 29(() ~ 9o\d£ = Y,9{ x i)Pj ~ £$(*i) a i ( 3 - 18 ) 

J ~ 1 j = 1 j=o 

N 2 

= - 9 (xo)a 0 + g(x N )a N - X)(~ 1 ) J »» _ ] 9{ x i) 

j=0 iV 1 

2 £. (-i YsM 


n 


where Cq = 2 = c^, Cj = 1 
so that 




(3.19) 


which proves the Lemma. 

However if one uses the new quadrature formula as in (3.5b) one gets 
Lemma (2.3): Let g Q be defined in (3.5a), then (3.9)(a) and (3.9)(b) are compatible. 
Proof: We start as in the last Lemma except now 

f {Pn- 29 ~ 9o)dt, — 9( x j)Pj ~ 9{ x j)ftj = 0 

J ~ 1 3 = 1 3 = 1 

which yields the result. 

We have thus established that the use of the Clenshaw- Curtis formula may lead to an 
approximation that does not satisfy (3.1b). The new quadrature formula alleviates this 
problem. 
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IV. NUMERICAL RESULTS 

In Section II we showed that the use of the Clenshaw-Curtis formula (2. 8), (2. 9) may cause 
problems when one attempts to solve the steady state problem corresponding to (2.1). In this 
Section we report on the numerical solution of equation (2.1), (2.2) with two different source 
functions, g(x ) = cos37rx and g{x) = cos57rx. To advance in time we use the fourth order 
Runge-Kutta scheme; a grid of 19 points was used. In table I we summarize the results for 
the first case: g(x ) = cos37rx. The first column gives the L 3 deviation from the steady state 
solution when (2.6a,b)-(2.7) is used, i.e. no modification for the right hand side. A linear 
growth is observed in the deviation from the steady state. The solution did not converge 
even after 200,000 time steps. 

In the second column we give the results for (2.6c), in which the approximation to the 
mean of g is obtained with the use of the Clenshaw-Curtiss formula. It is interesting to note 
that the results are essentially the same; the subtraction did not improve the situation. In 
the third column we present the results for (2.22a); convergence was obtained after 10,000 
time steps. 

In table II we present the same results for g{x) = cos5ttx. Basically, the results are the 
same as in table I; however, the divergence is much more rapid than in the first case. The 
new method converged after 15,000 time steps. 
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Table I 

L 2 Errors for g(x ) = cos 3 ttx N = 19 


No. Of 
Time Steps 

No 

Subtraction 

Clenshaw- Curtiss 
Subtraction 

GS 

Subtraction 

0 


IfiHI 

1.778 (-31 

10000 


IMBM 

2.669 (-7) l 


wmMSSl 


2.669 (-7)* 1 


1.492 (-5) 

1.494 (-5) 

same | 

150000 

2.234 (-5) 


same = 

200000 


III IMBB 

same 1 

Growth 

7.4165 (-7) 

7.4277 (-7) 

0 E 

■=- 


* converges after 25000 steps 


Table II 

Li Errors for g{x) = cos 5 ttx N = 19 


No. Of 
Time Steps 

No 

Subtraction 

Clenshaw- Curtiss 
Subtraction 

GS = 

Subtraction _ 

0 " 


6.320 (-4) 

6.320 (-4) 


1.4509 (-3) 

B— 

2.4405 (-4) E 

15000 

UftttMft)] 

■ III 1 MIIM 

2.4554 (-4) 

100000 



same m ._ 

150000 

■Wffl 

— M 

same — 

200000 

2.581 (-2) 


same 

Growth 

6.46 (-4) 

6.75 (-4) 

0 
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Finally, we show the results of a two-dimensional calculation for; 

u xx + Uyy = g{x,y) = cos -KX ■ cos iry 

with 

u x — 0 on (±l,y) (4T&) 

u y = 0 on (x,±l). (4-lc) 

Solution is carried out by the influence matrix technique, with the related Dirichlet Poisson 
problems solved via the direct solution method of tensor-product diagonalization. In order 
to set the level of the solution, that is, to set the arbitrary constant in the solution of the 
pure Neumann problem, one must set the value of the solution at one point on the boundary 
by replacing one row of the influence matrix. However, if the source function does not satisfy 
the discrete version of the compatibility relation 

J 1 J^g(x i y)dxdy = 0 (4.2) 

then large distortions appear in the solution near that point. On the other hand, if (4.2) is 
discretely satisfied, then the solution near the set point will be regular. 

In Figure 1 is shown isolines of the solution to (4.1) , for which g(x,y) was replaced by 

g(x , y) - £ £ 3( x h yk)ot k ai ( 4 -3) 

1=0 fc= 0 

that is, using the Clenshaw- Curtiss formula to satisfy compatibility. The location of the 
boundary point which is set in the influence matrix is obvious; the true solution is swamped 
by the artifact distortions. In Figure 2 is shown the isolines for the solution wherein the 
compatibility relation is satisfied using the quadrature formula developed here. The solution 
is now smooth, symmetric, and regular near the set point. 
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Fig. 2. Isolines of solution to (4.1) with GS subtraction. 
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